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Abstract: 

We present details of a lattice study of infrared behaviour in SU(3) gauge theory with twelve massless fermions in 
the fundamental representation. Using the step-sealing method, we compute the coupling constant in this theory 
over a large range of scale. The renormalisation scheme in this work is defined by the ratio of Polyakov loops in 
the directions with different boundary conditions. We closely examine systematic effects, and find that they are 
dominated by errors arising from the continuum extrapolation. Our investigation suggests that SU(3) gauge theory 
with twelve fiavours contains an infrared fixed point. 
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I. INTRODUCTION 



The origin of electroweak (EW) symmetry breaking is one of the most important research topics in physics. With 
the progress of experiments at the Large Hadron Colhder (LHC), it is urgent for a theoretical understanding for 
the mechanism of the mass generation and its relation to EW symmetry breaking. One appealing scenario for this 
mechanism is the technicolour models [H, ^ . These models involve new asymptotically-free gauge theories in which 
the coupling constants become strong at the TeV scale. The strong coupling can induce condensates to generate 
mass gaps, and asymptotic freedom leads to the absence of the naturalness problem. In order to dynamically 
suppress the flavour-changing neutral currents (FCNC), and to evade the constraints from precision EW data, 
it is important that the candidate theories exhibit the "walking" (quasi-conformal) behaviour and contain large 
anomalous dimension for the technifermion mass term . 

In recent years, there has been a significant amount of work in search of gauge theories viable for walking- 
technicolour model building. The most important task in this endeavour is the determination of the critical 
number of massless fermions, given the gauge group and the fermion representation, above which a theory is con- 
formal in the infrared (IR). For theories involving fermions in the fundamental representation, this is denoted as 
the critical number of flavours, N" . For N^'^ < Nf < (-^/^ is the number of flavours above which asymptotic 
freedom is lost), the theory contains an infrared fixed point (IRFP). A candidate walking-technicolour theory with 
fundamental fermions is believed to have the number of fiavours just below Nj". This makes the determination 
of iV™ a task with phenomenological significance, in addition to its importance in field-theoretic studies. Since 
the couplings must be strong at low energies in these theories, nonperturbative methods, such as the Schwinger- 
Dyson equation and gauge-gravity duality, have to be employed. Amongst these, lattice gauge theory is the only 
first-principle tool, and has been applied by many groups in this research avenue 

Of all the theories which have been investigated using the lattice technique, the value of A^" for SU(3) gauge theories 
with fundamental- representation fermions remains a controversy. Although several groups Irl- llli [26l. [33l |46| found 
evidence that SU(3) gauge theory with Nf = 12 is conformal in the IR, authors of Refs. [33, [S^l argued that 
chiral symmetry is broken in this theory. In this paper, we report our study of this theory, using the step-scaling 
method to compute the running coupling constant. We adopt the Twisted Polyakov Loop (TPL) scheme [49l - [5l| . 
This article complements the letter Q which was released in 2011 with other colleagues on this collaboration, and 
contains more details of our simulations and improved analysis using more data. In Ref. @ , we concentrated on 
the analysis with the step size, s, set to 1.5, while here wc emphasise the case in which s equals two. Furthermore, 
in the improved analysis with new data, as presented in this paper, we significantly reduce the correlation between 
data for the step-scaling functions on different lattice volumes. This makes the continuum extrapolation simpler 
and better controlled, compared to the analysis published in Ref. Q. We will discuss this in detail in Sec. IV CI 
Related to this work and Ref. Q, we have also published conference proceedings (s^ - ls^ . as well as for a similar 
project on SU(2) gauge theory with eight flavours (55| . 

In addition to computing the running coupling constant, we also obtain the ratio between the step-scaling function 
and the coupling constant, which becomes one when the /3— function is zero. To claim the discovery of the IRFP in 
an asymptotically-free gauge theory, we have to demonstrate that this ratio is indeed one in the ultraviolet (UV) 
and the IR, while being obviously different from this value between these two regimes. Our study suggests that 
SU(3) gauge theory with Nf = 12 contains an IRFP around the TPL-scheme coupling constant, 

gl - 2.0. (1) 

Amongst systematic effects that we estimate, errors arising from the continuum extrapolation dominate. We also 
notice that some of our procedures in performing this extrapolation lead to weaker evidence for the existence of 
the IRFP. Details of the estimation of systematic errors will be presented in Sees. fVl and I VII 

Our finding for the evidence of the existence of the IRFP agrees with the result of Refs. 0,lll, where the Schrodingcr- 
functional (SF) scheme [H,!!^] was used in defining the coupling constant, and the calculation was performed using 
the same gauge and fermion actions. The values of are different because of scheme dependence. Here we also 



^ There have also been many works on walking-technicolour model building using the gauge/gravity duality, as reviewed in Ref. [5 
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stress that the lattice artefacts can be very different in these two schemes. In particular, the SF scheme contains 
0{a) (a is the lattice spacing) lattice artefacts through the introduction of the boundary terms^, while in the TPL 
scheme the lattice artefacts remain of O(a^), making the continuum extrapolation more reliable. 

This paper is organised in the following way. In Sec. |lll we review twisted boundary conditions and the Twisted 
Polyakov Loop scheme. Section IIIII contains the details of our simulation strategy and parameters. We describe 
our analysis procedure in Sec. El give our results and discussion in Sec. IVIl and conclude in Sec. IVIII Appendix [K\ 
contains the study of the eigenvalue spectrum of the Dirac operator used in this work. Values of plaquctte and the 
raw data for the TPL-scheme coupling constants are presented in App. [B] 



II. TWISTED POLYAKOV LOOP SCHEME 



In this section, we give the details of our definition of the renormalised coupling constant in the twisted- Polyakov- 
loop (TPL) scheme [sol . [Slf . This scheme makes use of twisted boundary condition (TBC) [s^, which is implemented 
on the link variables, Ufj,(n) (/i = x, y, z, t is the Lorentz index and n is the position of a lattice site), through 

C/^(n + i)L,/a) = n,U^{h)ni, (2) 

where is the (dimensionful) box size in the ly direction (with j> denoting the unit vector), and a is the lattice 
spacing. The "twisting matrices", il^, act in the colour space. In this work, we apply TBC for ly = x,y, while 
maintaining periodic boundary condition (PBC) for the other two directions. This means 

n.^fit^ 1. (3) 

For SU(3), the twisting matrices, flx.y, satisfy 

n^nl = 1, (n^f = l, Tr [n^] = O, for ^l^x,y. (4) 
In this work, we explicitly implement [59| 

e-2"/3 

e2-^/3 I . (5) 
1 
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The inclusion of fermions is not straightforward when TBC, Eq. ©, is imposed on gauge fields. In order to 
maintain gauge invariance and single- valuedness of the fermion field, 'ijj{n + xL^/a + yLy/a), under the application 
of two boundary twistings (in x and y directions) with different orderings, it is necessary to introduce the "smell" 
degrees of freedom [60| . This quantum number is carried by fermions. The number of smells, Ng, is equal to the 
number of colours, Nc- Twisted boundary condition on fermion fields is given by, 

C(« + i^Lja) = e"/3r!eV^(n) (r!,)^„ , (6) 

where a and b are colour indices, and the twisting matrices, £7^, have been generalised to act on the smell degrees 
of freedom (indices a and 13). The factor e"/"^ is introduced only for v = x,y, to remove the zero-momentum 
modes in these directions. For v = z,t directions, we implement ordinary PBC, '0(n -I- vLnj a) = ijj{n). Since the 
smell quantum number is not carried by the gauge fields, it can be considered as additional flavours. Therefore the 
number of flavours in simulations involving dynamical fermions with TBC has to be a multiple of Ng {=Nc). 



^ In Rcfs. it was found that such 0(a) lattice artefacts can be numerically very small. 
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The Polyakov loops in the twisted directions, u = i, y, are 

W Ux {fix = j, fly, fiz, fit) 



Pxinx,fiy,fit) = Tr 
Py{fix,fiz-,fit) = Tr 



WUy {fix, fly = j,flz,flt) 



^ gi27rn^a/(3L^) 



(7) 



The extra factors outside the square brackets are introduced to maintain gauge and translation invariance. The 
renormalised coupling constant can be defined via the ratio between correlators of Polyakov loops in the twisted 
and periodic directions, 



(Px{nt=0)^Px{fit=Lt/{2a))) 
{P,{fit=0yp,{fit = Lt/{2a))) 



(8) 



where Pz^t are ordinary Polyakov loops in the directions with PBC. In this study, we always use hypercubic lattice 
Lx/a, = Ly/a = Lz/a = Lt/a = L/a. The proportionality factor k can be extracted by computing the above ratio 
in perturbation theory to 0{g'^), where op is the bare coupling constant. Using lattice perturbation theory, one 
obtains the lattice version of this factor [y] , 



(9) 



fc'^" = 0.03184 + 0.00453 +0 



The coupling, giatt, defined in Eq. ([8]) contains lattice artefacts, therefore depends on the lattice spacing as well as 
the volume. Its continuum-limit counterpart at fixed physical volume is defined as. 



gc = hm 5iatt, at fixed L. 



(10) 



The TPL scheme, as defined in Eq. ([8|), contains the feature that the renormalised coupling constant has the fixed 
value \/TJk ~ 5.6 in the IR limit {L oo). Therefore, in order to firmly establish the existence of the IR fixed 
point, we have to show that gc is significantly different from this value at the fixed point. 

Contrary to the SF scheme, the 0{a) lattice artefacts are absent in the TPL scheme. As explained in the following 
sections, it is important to control the continuum extrapolation in the step-scaling study of the running coupling 
constant. This makes the use of the TPL scheme very desirable. In Sec. |Vl we will show the lattice-spacing 
dependence of the TPL-scheme coupling constant. 



III. SIMULATION SETTING 



We give the details of our lattice simulation in this section. As discussed in Sec. [TTl the number of flavours in our 
calculation must be a multiple of Ng ~ Nc = 3. Since we are using staggered fermions, it also has to be proportional 
to the number of tastes, Nt = 4. In this work, we investigate the SU(3) gauge theory coupled to twelve flavours, 
which is allowed by these constraints. 



A. Step scaling 



Our goal is to measure the evolution of the running coupling constant over a wide range of scale. Given that the 
lattice imposes infrared (the volume) and ultraviolet (the lattice spacing) scales, the most convenient way to achieve 
this goal is the step-scaling technique. In this approach, we flrst measure the renormalised coupling constant, gutt, 
on the lattice in the scheme defined in Eq. ([5]). Since we perform computation at vanishing fermion mass, ^latt only 
depends on the lattice spacing and the lattice volume, L/a. Choosing a few values of L/a, we then simulate at a 



5 



wide range of /3 = &/g^, where go is the lattice bare couphng constant. This enables us to tune /3 (lattice spacing) 
to obtain the renormalised coupling in the continuum limit, 

5c {L) = giatt (/3l,£/ai) = 5iatt W2,L/a2) = . . . = glatt (/?noii/a«o) 7 (11) 

where rig is the number of choices oi L/a. Since (jc is independent of the lattice spacing, it is renormalised at the 
length scale L. In this work, we perform lattice simulations at 

L/a = 6, 8, 10. (12) 

Using the combinations of {(3, L/a) which lead to the same gc{L) (or u = g^), we compute the lattice step-scaling 
function, 

S(A,iM,w,s) = gfatt (A'SVa»)|„^g2^^j^^_i/„^) , (13) 

where i = 1, 2, . . . , np as in Eq. PT|) . and s is the step size. Since we can obtain uq results for S at the same physical 
volume, L, with different lattice spacings, this allows us to determine the continuum-limit step-scaling function, 

ct(w,s) = 5c 2(l) = lim S(/3,,i/ai,w,s). (14) 

V / a— 

In this work, we choose the step size s = 2, leading to the need for simulations performed on the lattice volumes, 

sL/a = 12, 16, 20, with s = 2. (15) 

For convenience, we define 

a{u) = a{u,s^2). (16) 

The step-scaling function is a scheme-dependent quantity, since it is simply the renormalised coupling constant 
computed at a certain scale. To facilitate a better method in demonstrating the existence of the IRFP, we compute 
the ratio 

rAu)^'-^. (17) 
u 

This ratio becomes one at the zeros of the /3-function. The existence of such zeros is independent of the rcnormal- 
isation scheme used in the calculation. In order to show that the gauge theory under investigation does contain an 
IR fixed point, we have to verify that r^riu) is one at both UV and IR regimes, while deviating from this value in 
between. 

A major source of systematic errors in the step-scaling method is the continuum extrapolation. It is a challenging 
task to properly address this issue. In order to have more information regarding this extrapolation and its possible 
systematic effects, we also perform simulation with L/a = 14, and resort to an interpolation procedure to obtain 
data for the TPL-scheme renormalised coupling on the lattice size L/a — 7. This enables us to carry out the 
investigation with 

{L/a = 6, 7, 8, 10) — > (2L/a = 12, 14, 16, 20) . (18) 

Here we stress that staggered fermions are used in this work, therefore it is not possible to have data directly on 
the L/a = 7 lattice. The interpolation procedure for obtaining such data is explained in detail in Sec. IV CI This 
interpolation in volume can introduce systematic effects, although it may result in more information regarding the 
continuum limit. Therefore, we only use the 4-point step-scaling analysis in Eq. psp as a means to estimate errors 
in the continuum extrapolation. 



B. Details of simulation parameters 



Our calculation is performed using the Wilson plaquette action for the gauge fields, and unimproved staggered 
fermions. We implement the standard Hybrid Monte Carlo (HMC) algorithm using the Omelyan integrator with 
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multi-time steps [6l|,[63 ■ To compute the inversion of the lattice fermion operator, biCGstab solver with convergence 
condition that the residue is smaller than 10^^® for molecular dynamics, and the accuracy of 10^^"* for the Metropolis 
tests, are used. A significant fraction of of our simulations were carried out on Graphics Processing Units (GPU's), 
where a mixed-precision solver with defect correction was implemented. The GPU codes were developed with 

CUBA m. 

To thermalise configurations in the Markov chains, we have used two procedures. In the first procedure, we start 
a simulation from a trivial gauge-field configuration 

Uf_,{nx,ny, h^,nt) = 1, (19) 

with fermion mass 

aruf - 0.5. (20) 

Then, we gradually decrease the mass to zero. In this process, we monitor the Polyakov loops in the untwisted 
directions, and make certain that the imaginary parts are non-vanishing. This ensures that the Markov chains 
progress mostly near the true vacua @. In the second procedure, we start with a configuration, 

U,in,,hy,n, = l,ht) = e-2W3 ^ Ut{n,,hy,h„nt = 1) - e+'''^^^ , (21) 
Uf_i{'hx,ny,nz,nt) ~ 1 elsewhere, 

which always results in non-zero imaginary parts in the Polyakov loops in the untwisted directions. It also produces 
the largest gap in the vicinity of zero in the fermion matrix. In this case, we can start the simulation directly with 
zero fermion mass, making this procedure significantly more efficient than the one implemented with the initial 
conditions of Eqs. (|19p and (|20p . In both cases, we observe that the simulations always stay near the true vacua, 
and tunnelling amongst these vacua occur occasionally. We will discuss this issue in more detail in Sec. IIVBI 

In order to implement the step-scaling investigation of the running coupling constant as discussed in Sec. IIII A[ we 
carry out simulations at the lattice volumes, 

L/a = 6, 8, 10,12,14,16,20. (22) 

For each volume, we simulate at several /3 values between 4 and 99, in the gauge action. Since the running is 
expected to be slow in SU(3) gauge theory with twelve fiavours, this large range of /3 is necessary to trace the 
coupling constant from the UV to the IR regimes. We aim at determining the Polyakov-loop correlators with 
statistical errors around 2.5% or smaller. For this purpose, a significant amount of gauge-field ensembles have to 
be generated. The raw data for the TPL-scheme renormalised coupling, as defined in Eq. ([5]), are given in App. [Bl 



IV. PLAQUETTE, POLYAKOV LOOP AND THE VACUUM STRUCTURE 

In this work, we perform several detailed checks on the simulations, in order to ensure that we are estimating 
the autocorrelation and performing the continuum extrapolations reliably. These checks include the lowest-lying 
eigenvalue spectrum of the Dirac operator, the plaquette values, and the phases of the Polyakov loops. The 
computation of the lowest-lying eigenvalues is presented in App. \^ while in this section we address the other two 
topics. 



A. Plaquette 

As shown in App. [Bl some of our simulations arc performed at small f3 values (coarse lattice spacings). It is 
necessary to check that these simulations are still in the weak-coupling phase, in order to make certain that at 
these /3 values, the theory is still in the same universality class as that with high—/? (fine lattice spacings). This is 
essential in order to ensure that the continuum limit can be reliably taken in our calculations. 
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FIG. 1: Upper panel: The values of plaquette as a function of "WCE" stands for "weak coupling expansion". Lower 

panel: The difference between the data points and the function, Eq (|24|l . fitted to the L/a — 20 data points only. 



For the above purpose, we examine the expectation values of the plaquette for many of our HMC simulations. The 
results are summarised in Table HII] in App. \S\ These expectation values are plotted in the Upper panel of Fig. [U 
where we also show the predictions from the weak coupling expansion for pure Yang-Mills theory, 

2 

plaquette ~ 1 ^ ^ (weak coupling expansion). (23) 

By comparing our data with this function, it is evident that all our simulations are in the weak-coupling phase, 
and are safe from being in the novel phase observed in Ref. [3^. We have also studied the volume dependence 
of the plaquette, by first fitting the data obtained on the largest lattice, L/a — 20, to a weak-coupling expansion 
formula (pi are the fit parameters), 

m-Po + j + ^ + ^ + ^,+^, (24) 

then computing the difference between the data points to this curve. The result of this investigation is shown in 
the lower panel of Fig. [1] This shows that finite-size effects are minor in the computation of the plaquette in this 
work. 



^ We thank David Schaich for private communications regarding this issue. 
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FIG. 2: Complex values (left panels) and the ratios between the real and imaginary parts (right panels) for Polyakov loops in 
the twisted (upper panels) and untwisted (lower panels) directions, in the first 25000 trajectories in the simulation performed 
at /3 = II. 15 and L/a = 16. 



B. Polyakov loops and vacuum tunneling 



The study of the plaquettes in the last section confirms that our simulations have been carried out in the weak 
coupling phase. In this phase, as pointed out in Ref. @, the true vacua in SU(3) gauge theory with fermions 
are always those in which the vacuum expectation values of the Polyakov loops in the untwisted directions are 
non- vanishing and complex. On the other hand, in the vicinity of the false vacua, the untwisted Polyakov loops are 
real. As for the Polyakov loops in the twisted directions, we expect that they will scatter around zero, configuration 
by configuration. 

Markov chains in our simulations can be trapped in the false vacua. However, by using the above property of the 
untwisted Polyakov loops, we can monitor the simulations and ensure that they arc mostly progressing near the 
true vacua. 

Investigating the Polyakov loops trajectory by trajectory, we first confirm that in the twisted directions, they are 
fluctuating around zero for all simulations in this work. This is shown for two typical cases in the upper panels of 
Figs. [Hand El 

Next, we study the Polyakov loops in the untwisted directions. In all our simulations, their values are non- vanishing 
and complex in all trajectories. The complex phase fluctuates around ±27r/3, indicating that the Markov chains 
are progressing near the true vacua. The lower panels of Fig. [2] demonstrate a case {L/a = 16, (3 = 11.15) in which 
the simulation stays near the vacuum with the phase of Polyakov loop being — 27r/3. For simulations performed 
at smaller L/a (fewer total degrees of freedom) and larger /3 (stronger coupling), tunnelling between the two true 
vacua may occur. One of such cases is shown in the lower panels of Fig. [3] Every time this takes place, we then 
investigate the Polyakov loop correlators trajectory by trajectory, ensuring that these correlators do not exhibit 
any "discontinuous" behaviour when the tunnelling happens. In Fig. |4l we show the result of this study for the 
corresponding simulation presented in Fig. |31 From these plots for the Polyakov loop correlators in the twisted 
and untwisted directions, we conclude that tunnelling between the true vacua does not result in artefacts which 
complicate the estimation of autocorrelation time. 
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FIG. 3: Complex values (left panel) and the ratios between the real and imaginary parts (right panel) for Polyakov loops in 
the twisted (upper panels) and untwisted (lower panels) directions, in the first 25000 trajectories in the simulation performed 
at /3 = 5.53 and L/a — 8. 
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FIG. 4: Polyakov loop correlators in the twisted (left panel) and untwisted (right panel) directions, in the first 25000 
trajectories in the simulation performed at /? = 5.53 and L/a = 8. 



V. ANALYSIS DETAILS 



In this section, we explain the details of our analysis. The statistical analysis in this work is performed using the 
bootstrap procedure, in which 1000 bootstrap samples are generated for each {L/a, 13). 



A. Autocorrelation and data binning 

As presented in App.|Bl we perform our calculations with a large number of HMC trajectories. The first step in our 
analysis is the binning of the raw data. In order to make certain that the binning procedure is reasonable, we study 
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the autocorrelation of the ratio, appearing in the left-hand side of Eq. (|H]), between the Polyakov loop correlators. 
To describe our investigation, we start from the autocorrelation function of primary quantities, [gS-IggI 



N-T 



r^M^) = ^ E i^^(') - (^^(^ + ^) - ^p) ■ (25) 

1=1 

Here, a and /3 label the types of primary quantities. In our case, ©1(1) and 02{i) arc Polyakov loop correlators 
of the i-th sample in the twisted and in the periodic directions, respectively. The quantity (5 5, is the average of 

By using F^^, the autocorrelation function of the Polyakov loop ratio, as in the left-hand side of Eq. ([8]), can be 
written as, 

2 

r(T) = E /"//sr^M^) (26) 



with, 



f ^^[9l\^± f ^^f9l]=^9l (27) 



We define the normalised autocorrelation function, 

r(0) 



Pi^) = -F77TT (28) 



which is normally assumed to behave as, 

p(r)-e"^. (29) 
The quantity ta is the autocorrelation time of single exponential autocorrelation. 

Since the integrated autocorrelation function is less noisy than p(r), wc use it to estimate the autocorrelation time 
between Polyakov-loop-corrclator ratios. Upon integrating over r, we obtain 



J p{T')dT' TA (^1 - e TA when r > ta. 



(30) 



The single-exponential form in Eq. (|29p is often a poor approximation to p(t), when the system contains degrees 
of freedom that are characterised by very different autocorrelation times. In general, the autocorrelation function 
can be multi-exponential, 

p(t) ^ E «fe ' ^ith E «fe = 1- (31) 

k k 

The integrated autocorrelation is 



r p{r')dr' ^Y.^'^^^ ^'^ (l-e-/^^ 
Jo I. ^ 



(32) 



k 

This function reaches a plateau t'^'^ au when t ^ t''^^ for all k. We use this criteria for the estimation of 

autocorrelation without explicitly determining t'^^'^ and ak- A more detailed study of autocorrelation times for 
conformal field theories will be reported in a separate paper [63|. 

In our numerical calculation, the integrated autocorrelation is defined as, 

©(^) = J + E P^^')- (33) 
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FIG. 5: Representative plots for the integrated autocorrelation of the ratio of Polyakov loop correlators at various values of 
L/a and /3. The physical volume increases from the top-left to the right-bottom corners. 



To estimate error in 9(r), we apply the Madras- Sokal formula [68l |. 

(Ae(r))2 - ^©M'- (34) 

Figure [5] shows Q{t) for the representative cases in this work. The separation between two decorrelated trajectories 
can be estimated by investigating the plateau of 6(t). As demonstrated in Fig.[5l this separation depends on the 
physical volume, L. It is around 20 on the smallest volumes, and about a few hundred to 1000 on the largest 
volumes. 

In App.IbI we show the details for the numbers of HMC trajectories in our simulations. For each choice of {L/a, (5), 
we divide the trajectories evenly into ^ 200 bins by averaging over them in each bin. These bins are then used to 
create 1000 bootstrap samples. From the result presented in this section, it is evident that our bin sizes are large 
enough compared to the autocorrelation times. This ensures that the data amongst these bins are decorrelated. 
We have also confirmed this with the Jackknife analysis using these and larger bin sizes. The statistical errors in 
this approach are almost the same as those in our bootstrap analysis, and they are stable against the change of 
the bin sizes. Figure |5] shows some examples for this Jackknife check for the TPL-schcme renormalised coupling 
computed at various /3 values on the L/a = 20 lattice. In this figure, each group of data points contains the results 
of the jackknife analysis with the numbers of bins set to 100, 200 and 500. The red point in the centre of each 
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FIG. 6: Dependence of the TPL-scheme renormalised coupling on the number of bins at various /3 values on the i/a = 20 
lattice. Each group of data points contains the results of the jackknife analysis with the numbers of bins set to 100, 200 and 
500. The red point in the centre of each group is the result of using 200 bins, as chosen in our bootstrap procedure. The /3 
values for the two blue points are slightly shifted for the purpose of presentation. They correspond to choosing the numbers 
of bins to be 100 (left) and 500 (right). 



group is the result of using 200 bins, as chosen in our bootstrap procedure. The /3 values for the two blue points 
are slightly shifted for the purpose of presentation. They correspond to choosing the numbers of bins to be 100 
(left) and 500 (right). From these plots, it is apparent that having 200 bins leads to enough trajectories in each 
bin, in order to correctly estimate statistical errors. 



B. Interpolation in /3 (bare coupling constant) 



In the step-scaling study of the running coupling constant, we first have to perform the tuning of the j3 values in 
Eq. (|lip . for the lattice volumes L/a = 6, 8, 10. In principle, this can be achieved by repeatedly adjusting (3 and 
carrying out new simulations, until Eq. (|lll) is satisfied to high accuracy. As discussed at the end of Sec. IIII Al we 
also want to obtain the TPL-scheme renormalised coupling on the L/a ~7 lattice through interpolation in volume, 
in order to estimate systematic errors in the continuum extrapolation. For this purpose, one has to tunc a different 
set of /3 values for L/a ~ 6, 8, 10 and interpolate to L/a = 7 at each step of this tuning. 

The above procedure is very time-consuming, and becomes impractical for studies in which one has to trace the 
coupling constant across a large range of length scale. This is the case in the current work. Therefore we resort to 
a variation of the above method. That is, we simulate at many /3 values for each L/a, and perform interpolations 
in /3 for the renormalised coupling constant, volume by volume. The choices of these /3 values are presented in 
App. [B] The use of this interpolation method inevitably introduces systematic effects in our calculation. We will 
address this issue in this section. 

Since we are simulating at a large range of bare coupling constant, it is a challenging task to have a well-inspired 
interpolation function in /?. One reasonable way to proceed is to note that in the large—/? (small bare-coupling) 
regime, one-loop perturbation theory has to be valid, and therefore at fixed L/a, 

"iatt = ffiitt(/3,i/«)«|=5o (for/?»l), (35) 

where is the bare gauge coupling. This motivates the use of polynomial functions in 1//3 to perform the 
interpolation. Since we have data for many /3 values (see App. [B| for each L/a, it is in principle possible to have 
high degrees of polynomials for these fits. Such high-degree polynomials will generally fit all the data points. On 
the other hand, the Runge phenomenon may occur in this procedure, resulting in artificial oscillatory behaviour 
of the fit functions. In order to avoid this artefact in the /3— interpolation, we note that the renormalised coupling 
should always be non-decreasing with growing lattice spacing (i.e., decreasing p) at fixed L/a, otherwise the theory 
will be in the strong-coupling phase and the continuum extrapolation cannot be reliably performed. 
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From our study of the plaquettes in Sec. IIV Al it is evident that our simulations are all carried out in the weak- 
coupling regime. This is also reflected on the data points plotted in Fig. [71 in which we see that all our renormalised 
couplings, uiatt, arc non-decreasing when /? decreases. This leads to the use of the non-decreasing polynomial, 



/(«o) = j dua c,„u^" = (^li^^e uo=- = ^-^y (36) 



"latt 



in the /3— interpolation procedure at fixed L/a. We implement the constraint from perturbation theory, Eq. ([55)1 . 
which results in, 

/lo = 0, /ii = 6 (then cq = \/6). (37) 
This constraint leads to the number of fit parameters, 

^param = A^dcg = — , (38) 

where iVdog and Nh are defined in Eq. (|36p . The use of the non-decreasing polynomial ansatz makes the Runge 
phenomenon milder compared to the simple polynomial fits. The inverse of the fit function in Eq. (|36p is also 
single-valued. This is essential in the step-scaling method. The results of applying this (uncorrelated) fitting 
procedure in the /3 interpolation are shown in Fig. [T] The optimal choices of A^parami leading to the best (smallest) 
X^/d.o.f., are listed in TablelH 

In order to estimate systematic error resulting from the interpolation in /3, we change the fit function from Eq. p6p 
to a simple polynomial function, 

Wdcg 

Wlatt = /(wo) = ^ CmU™, (39) 
m=0 
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Non-decreasing Polynomial 


L/a 


^ par am 


Y^/d o f 


6 


7 


1.654259 


8 


5 


0.837240 


10 


5 


0.828201 


12 


4 


1.597743 


14 


4 


2.498352 


16 


4 


0.834323 


20 


7 


0.685983 



TABLE I: Left: The x^/d.o.f. of the /3 interpolation using Eq. 
Right: The x^/d-o.f. of the /3 interpolation using Eq. ((39|. iVp 



Simple Polynomial 



L/a 


-/Vparam 


, 2 / 1 „ r 

X /d.o.l. 


6 


8 


1.580600 


8 


11 


0.652351 


10 


5 


0.819650 


12 


4 


1.612676 


14 


6 


2.608492 


16 


4 


0.837765 


20 


7 


0.689820 



'^'^„ ^ is the number of fit parameters. 



~ Ndcg, — 1 is the number of fit parameters. 



with the constraint, 

£o = 0, Si = 6. (40) 
from the vahdity of perturbation theory at high—/?. This constraint results in the number of fit parameters, 

iVparam - ^dcg " 1- (41) 

The values of A^param for the best x^/d.o.f. are presented in Table ID 



C. Interpolation for L/a = 7 



As indicated at the end of Sec. IIII A[ it is desirable to gain more information regarding systematic errors in the 
continuum extrapolation. In view of the fact that our reference input renormalised couplings are computed on 
L/a = 6, 8, 10, a practical way to proceed is to have data for L/a = 7. This enables us to attempt the step-scaling 
study, 

{L/a = 6, 7, 8, 10) — > [sL/a = 12, 14, 16, 20), where s = 2, (42) 
without having to perform simulations on large lattices, such as L/a = 24. 

Since staggered fermions are used in this work, we have to use an interpolating procedure to obtain uiatt for 
L/a = 7. To have a well- motivated method for this interpolation, we resort to the /?— function of the theory. It 
is well-established that the coupling constant in SU(3) gauge theory with twelve flavours runs slowly compared 
to, e.g., QCD. This is reflected on the fact that a small change in the renormalised coupling has to result from a 
significant variation of the scale. As shown in Fig. [71 this is indeed the case. Namely, enlarging the box size by 
a factor of two induces very little changes in uiatt, and one can locally approximate the /3— function using a linear 
form 

^rfuiatt ^ fj{u^^^^.) Kiai+bi Uiatt, (43) 

dL 

where ai and bi are unknown parameters. We stress that this approximated form is not based on perturbation 
theory, and is only valid within a small range of wiatt- That is, in different ranges of wiatt, the parameters, ai and 
6/, have different values. 

To determine uiatt on the L/a = 7 lattice, we use our data on the L/a = Q, 8, 10, 12 lattices, and interpolate with 
the function, 

wiatt = Al+Cl (-\ , (44) 
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at fixed lattice spacing. The unknown coefficients, Al, Bl, and Cl are related to ai and bi, and the integration 
constant in solving Eq. p3|) . Figure [8] shows two representative plots for the interpolation using Eq. p4|) . It is 
obvious that the interpolation is smooth, and the values of the coefficients, A^, and Cl, can vary significantly 
in different ranges of itiatt- The fits presented in Fig. |S] are performed on the L/a = 6,8,10,12 data without 
^—interpolation. 

Equations (|43)) and (|44)) are used to motivate an interpolation function in i/a at fixed a {f3 value). However, the 
the effects of the lattice spacing can appear as powers of (a/L)"^ in our simulations. This means the data points used 
in each of this volume interpolation may have different lattice artefacts, leading to systematic effects introduced 
in this procedure. In view of this, we do not include the L/a = 7 data in our central analysis procedure, and 
only use them to perform the step-scaling investigation in Eq. (|42p as a means to estimate errors in the continuum 
extrapolation. 

Another issue in this volume- interpolation method for obtaining the L/a ~ 7 data is statistical correlation. The 
procedure is carried out using uncorrelated fits in this work. However, it is natural to expect that there will be 
correlation between the L/a = 7 (interpolated) data and those extracted directly from independent simulations on 
L/a = 6, 8, 10, 12. This correlation has to be closely examined, since all these data are used in the investigation of 
the continuum extrapolation, as discussed in Sec. IV Dl For this purpose, we study the likelihood function, 

L{lH, Uj) = , CXp \^{U^- Ui) [COV^^] {Uj -Uj) \ , (45) 

ZTTy'det (Gov) L ^ J 

where Ui denotes uiatt computed on the lattice volume L/a = i. Here is kept as a variable, and Ui is its central 
value for this quantity from our simulation. The symbol Gov is the covariancc matrix which can be computed from 
the bootstrap samples of Ui and Uj obtained from numerical calculations. 

Our investigation shows that, although the coupling constant on the L/a ~ 7 lattice is interpolated using those 
on the L/a ^ 6,8, 10, 12 lattices, it only shows significant correlation with that on the L/a ^ 8 lattice. In Fig. ^ 
wc display an example of this likelihood- function study performed for /3 — 5.53. It is obvious from this figure that 
coupling constants on L/a = 7 and L/a ~ 6, 10, 12 exhibit very small correlation, while it is the opposite between 
L/a = 7 and L/a = 8. The corresponding covariancc matrices for the example in Fig. [3] are 
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FIG. 9; The likelihood function plotted against (ug, mt), (itg, u-j), (mio, uy), and (1112, U7) at /3 = 5.53. The dashed curves 
indicate the standard error ellipses. 



000235 





000104 


000104 





000999 


005731 





001473 


001473 





000999 


005945 
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000776 





000999 
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000276 





000999 




Coy(^=5.53) ^ 

y 0.000776 0.000999 j 

Cov(^- -) ^ / 0.005090 0.000276 \ 
\ 0.000276 0.000999 ) 

The volume, L/a = 7, is one of the "small" lattices, on which we compute the reference coupling instead of the 
step-scaling function. The importance of the above study is the demonstration that there is negligible correlation 
between data on this lattice and that on the "large" lattice, L/a = 12, from which we compute the step-scaling 
function. Furthermore, the statistical errors of the data obtained on all our small lattices are small. In view of 
this, it is reasonable to expect that this correlation between the L/a = 7 and L/a = 8 TPL coupling constants 
does not necessitate correlated fits in the continuum extrapolation. 
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The above study of the data correlation also leads to the conclusion that one has to be very cautious about 
interpolating wiatt in L/a. In certain analysis procedures, such as the one we adopted in Ref. @ by setting the step 
size to 1.5, large correlation amongst data used in the continuum extrapolation can occur. 



D. Continuum extrapolation for tlie step-scaling function 

The last step in our analysis is the continuum extrapolation for the step-scaling function, defined in Eqs. ()14p 

and (|16p . Since unimproved staggered fermions and the Wilson plaqucttc action are used in this work, we will 
investigate (a/L)^ dependence in the lattice step-scaling function, Y,{(3,L/a,u,s = 2). In Fig. [101 this dependence 
is displayed at representative values of u in the regimes of weak, intermediate and strong coupling. From this figure, 
it is obvious that effects of the lattice artefacts grow with increasing u, as expected. In the region u < 0.8, we see 
that the step-scaling functions show insignificant dependence on the lattice spacing, and are almost consistent with 
the input reference coupling. On the other hand, in the strong-coupling regime, the a— dependence in S becomes 
noticeable, necessitating good control of the continuum extrapolation in the investigation of the existence of the 
IRFP. It is worth noting that the lattice artefacts tend to make the step-scaling function larger than its continuum- 
limit counterpart, especially in the strong-coupling regime. This feature is different from what was discovered in 
the Schodinger-functional scheme 0, § . 

In performing the continuum extrapolation for our central analysis procedure, we use our simulation results for 
E(/3, L/a, u, s = 2), obtained at sL/a = 12, 16, 20, and carry out the linear fit {(Ti{u) and Ai are the fit parameters), 

E(/3, L/a, u, s = 2) = ai{u) + Ai{u) (|-)' , (47) 

with the /3— values for various L/a determined by tuning the coupling, u, to be the same on the corresponding 
small lattices {L/a = 6,8,10). This procedure does not include the L/a = 7 data which are extracted with an 
additional volume-interpolation, as detailed in Sec. IV CI 

To estimate systematic errors in the continuum extrapolation, we include the volume-interpolated, L/ a ~ 1 data, 
as well as the step-scaling functions computed on the lattice, sL/a = 14. We first perform the quadratic fit {aq{u), 
Ag and Bq are the fit parameters), 

E(/3, L/a, u, ,s = 2) = a„iu) + Aq{u) (^)' + Bq{u) (^)' , (48) 
to implement the 4-point step-scaling method in Eq. (|18p . 

In order to further account for systematic effects arising from the continuum extrapolation, we perform two addi- 
tional linear fits: 

1. Using the data for the step-scaling functions from sL/a = 14, 16, 20 {L/a = 7, 8, 10). 

2. Using the data for the step-scaling functions from sL/a ~ 12, 14, 16, 20 {L/a — 6, 7, 8, 10). 

Figure [TT] shows representative plots of the continuum extrapolation using the above procedures (quadratic fit and 
the three linear fits). From these plots, we observe that a; and aq are well consistent with each other at intermediate 
and strong couplings. In the weak-coupling regime (top row of Figure [TT]), we notice that the quadratic fit, and the 
3-point linear fit using the i/a = 7, 8, 10 data are not consistent with the other two procedures. They result in a{u) 
smaller than u after the continuum extrapolation. However, we stress that in this regime, the lattice step-scaling 
function, E, demonstrates very mild lattice-spacing dependence, and is almost consistent with the input reference u. 
This is the consequence of asymptotic freedom. Furthermore, our data do not show significant 0{a^) contributions 
in the continuum extrapolation at strong and intermediate couplings (center and bottom rows of Figure fTTjl. where 
the lattice artefacts are expected to be larger compared to the small— u region. In view of this, we conclude that 
the quadratic fit, and the 3-point linear fit using the L/a = 7,8, 10 data can be artificially amplifying statistical 
fluctuations and leading to unreliable results in the weak-coupling regime. In order to properly address this issue, 
one has to generate data with very high statistical accuracy {e.g., < 0.5%) at large /3 values. This is beyond the 
scope of this work, since our main focus is on the existence of the IRFP in the strong-coupling regime. 
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FIG. 10: Lattice-spacing dependence of the step-scaling function (SSF) in weak, intermediate and strong coupling regimes 
(from the top). The horizontal lines indicate the central values of the input reference u. 
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FIG. 11: Representative cases of the continuum extrapolation for the step-scaling functions using the procedures discussed 
in the main text. The 3-point linear extrapolation using data on L/a = 6, 8, 10 is the central procedure. The horizontal lines 
indicate the central values of the input reference u. As discussed in the main text, the quadratic fit, and the 3-point linear 
fit using the L/a = 7, 8, 10 data can lead to unreliable results in the continuum limit in the weak-coupling regime (top row). 



VI. FINAL RESULTS AND DISCUSSION 



In this section, we present the final results of our analysis, and discuss the estimation of systematic errors. We 
begin by showing the result from our central analysis procedure for the ratio r„ = a{u)/u, defined in Eq. (jl7p . 
In performing this central-procedure analysis, we first interpolate in the bare coupling, /3, for simulation data 
obtained at L/a = 6, 8, 10, 12, 16, 20, using the non-decreasing polynomial function in Eq. ([55)1 with the constraint 
from Eq. ([57]) . and the polynomial degrees and the numbers of fit parameters presented in Table HI We then carry 
out the step-scaling of 

L/a = (6, 8, 10) — > 2L/a = (12, 16, 20), (49) 

by extrapolating the step-scaling function to the continuum limit with the linear form in (a/L)^, Eq. (1471) . Result 
of this central analysis is shown in Fig. 1121 which demonstrates evidence for the existence of an IRFP. 
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FIG. 12: raiu) from the central procedure. 



1.2 
1.15 

1.1 
1.05 
1 

0.95 
0.9 
0.85 



1.2 
1.15 

1.1 
1.05 

1 

0.95 
0.9 
0.85 



2-loop 



0.5 1 1.5 2 

Reference Coupling u 



2.5 



2-loop 



1 1.2 1.4 1.6 1.8 2 2.2 2.4 
Reference Coupling u 



1.2 
1.15 

1.1 
1.05 
1 

0.95 
0.9 
0.85 



2-loop 



1 1.2 1.4 1.6 1.8 2 2.2 2.4 
Reference Coupling u 




1 1.5 2 
Reference Coupling u 



FIG. 13: Plots for ra{u) obtained from our procedures for estimating systematic errors. Top left: r^iu) from simple 
polynomial interpolation in /3, Eq. (|39|) . Top right: rcr{u) by performing the continuum extrapolation using quadratic 
function in (a/L)^. The rest is the same as the central procedure. Bottom left: raiu) by performing the continuum 
extrapolation using linear function in (a/L)^, with L/a — 7,8,10. Bottom right: ra-{u) by performing the continuum 
extrapolation using linear function in (a/L)'^, with L/a = 6, 7, 8, 10. 



Next, we discuss the estimation of systematic effects arising from the /3— value (bare-coupling) interpolation and 
the continuum extrapolation. For this purpose, we perform the changes in the central procedure. These changes 
are carried out independently, i.e., we vary one component in the central procedure, while keeping the other fixed. 

We begin by varying the interpolation in the central procedure. This is carried out by changing the non- 
decreasing fit function in Eq. ([55]) . to the simple polynomial form in Eq. ([M)) with the constraint of Eq. PO]) and 
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the numbers of parameters reported in Table HI The result of this procedure is shown in the top- left panel of Fig. 1131 

In order to estimate systematic errors associate with the continuum extrapolation, we perform various fits with the 
inclusion of the L/a = 7 and L/a = 14 data. First, we perform the quadratic fit using Eq. (pS]). This leads to the 
result for r„ (u) as depicted in the top-right panel of Fig. 1131 As expected, this extrapolation strategy results in 
large statistical errors. In addition to the quadratic fit, we also carry out the two linear continuum extrapolations 
discussed in Sec. IV D[ 

L/a =(7, 8, 10) — > 2L/a= (14, 16,20), (50) 
L/a = (6,7,8, 10) — > 2L/a = (12, 14, 16, 20). 

The result from the first these procedures is presented in the bottom-left panel of Fig. [T31 while that from the 
second one is shown in the bottom-right panel of Fig. [T^l As discussed at the end of Sec. IV Dl the quadratic 
fit, and the 3-point linear fit using the L/a ~ 7,8,10 data can lead to unreliable continuum extrapolations in 
the weak-coupling regime. Therefore, in Fig. [13] we only show results at intermediate and large u for these two 
procedures. 

In Fig. [T^l and in the top-left and the bottom-right plots in Fig. [TSl it is observed that these procedures lead 
to rcr(u) consistent with one in the UV and the IR, while statistically different from this value between these 
two regimes. This suggests that there exists an IRFP in SU(3) gauge theory with twelve flavours. However, the 
continuum extrapolations using 4-point quadratic fit (the top-right plot in Fig. I13p and 3-point linear fit without 
the L/a = 6 data (the bottom-left plot of Fig. [T5)) lead to weaker evidence for the IR conformal behaviour. For 
these two procedures, in addition to the difficulty in the continuum extrapolations in the weak-coupling regime 
(discussed at the end of Sec. IVPp . we also observe large errors in the IR regime, leading to no apparent feature 
that r^{u) crosses one. This phenomenon is actually the consequence of the "double crossing" behaviour in some 
bootstrap samples. Namely, in these samples, r^{u) crosses one from above, and then turns around to cross the 
same value from below in a slightly larger u. We stress that out of 1000 bootstrap samples we have created in 
this work, r^riu) in more than 680 (Icr) of them cross the unity from above in the IR regime, when the continuum 
extrapolations are performed with the quadratic fit or the 3-point linear fit without the L/a = 6 data. This leads 
to hints of the existence of an IRFP using these analysis procedures. In order to illustrate this point, in Fig. [HI we 
plot 100 bootstrap samples in the intermediate— and strong— u regions in these two procedures. 

In Table[TTl we summarise the values of obtained from the above procedures. Because it is challenging to precisely 
estimate systematic effects, as discussed above, we take a conservative approach to conclude that in SU(3) gauge 
theory with twelve flavours, our data suggest the existence of an IRFP around, 

5*' - 2.0. (51) 

This result is similar to what we obtained with other collaborators using a different analysis procedure @ by setting 
the step size to be 1.5. Here we stress that it is more challenging to control systematic effects and the correlation 
amongst data points in the procedure in Ref. , because of the need for many interpolations in lattice volumes 
when computing the lattice step-scaling function, S. 
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P interpolation 


continuum extrapolation 




non-decreasing polynomial 


3 point linear, L/a — 6, 8, 10 


2.02(18) 


simple Polynomial 


3 point linear, L/a — 6,8,10 


2.02(21) 


non-decreasing polynomial 


4 point linear 


2.06(15) 


non-decreasing polynomial 


3 point linear, L/a = 7, 8, 10 


> 1.66 


non-decreasing polynomial 


4 point quadratic 


> 1.62 



TABLE II: g'i from various procedures. The first row describes the central procedure. 

The result in Eq. (|51|) is much smaller than that obtained in the SF scheme 0, ^ , 

(.9r))'-4.5. (52) 

The significant difference clearly indicates that the two schemes are very different. It should also be noted that in 
the TPL scheme, there is an upper bound for the renormalised coupling constant, as discussed in Sec.|TTl This may 
result in slower running behaviour compared to the SF scheme. 



VII. CONCLUSION 



In this paper, we present our work on the lattice study of IR behaviour in SU(3) gauge theory with twelve flavours. 
We use the step-scaling method to investigate the running coupling constant over a large range of scale. Our 
renormalisation scheme is defined via the ratio of Polyakov loop correlators in the twisted and untwisted directions. 
In particular, we compute the ratio, r^iu) defined in Eq. (|17|) . between the step-scaling function and the input 
renormalised coupling. In our central analysis procedure, we perform the continuum extrapolation using the 3-point 
linear fit with the reference coupling computed on the L/a = 6,8,10 lattices. Data on these lattices are free of 
volume interpolation. Using this procedure, we find that this theory contains an IRFP at around gl ^ 2. 

In this work, we have investigated systematic errors in the bare coupling interpolation for each lattice volume, and 
the continuum extrapolation. We have performed reasonable variations on these interpolation and extrapolation, 
and carefully examined possible correlation amongst data points used in the continuum extrapolation. We find that 
the dominant systematic effect arises from the continuum extrapolation. To gain information about possible errors 
in this extrapolation, we compute the step-scaling function on the L/a = 14 lattice, obtain the reference input 
renormalised coupling for the L/a = 1 lattice using an interpolation procedure, and then study the continuum 
limit using the 4-point linear and quadratic fits, as well as the 3-point linear fit without the L/a = 6 data. We find 
that all our analysis procedures result in evidence for the existence of an IRFP, although the latter two continuum- 
extrapolation methods result in significant errors. In view of this, the result of our work suggests that SU(3) gauge 
theory with twelve fermions in the fundamental representation contains an IRFP. 

Our finding shows that the conformal window for SU(3) gauge theories with fundamental fermions rr iay lie below 
N j = 12. Although this conclusion agrees with most other studies 04111 HI, IH, , the result in Rcf. [Ijjll^l leads 
to the opposite conclusion. Combining this information with the recent result from the Nf = 10 calculation [T2j . 
this can indicate that the N f = 12 is already very close to the lower bound of the conformal window. 
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FIG. 15: Plots of the 16 lowest-lying positive eigenvalues (in lattice units) of the staggered fermion operator in this work. 
The two plots show a coarse (left) and a fine (right) lattices of the same physical volume. Each line in these plots connects 
the eigenvalues computed on the same gauge configuration. 
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Appendix A: Low-lying eigenvalues of the Dirac operator 



To check the effects of taste-symmetry breaking in staggered fermions, we study positive low-lying eigenvalues of 
the Dirac operator. In this section, we present a typical case of taste-symmetry restoration when approaching the 
continuum limit at fixed physical volume. For this purpose, we compare the following two cases: 

1. L/a =10, 13 = 20.13, in which gf^^^ = 0.4031(76). 

2. L/a = 20, /3 = 20.00, in which gf^^^ = 0.4064(76). 

The renormalised coupling for these two cases, as shown in the table for the raw data in App. |Bl are well consistent 
within statistical error. This means that the physical volumes are almost the same, while the lattice spacing of the 
second case is half of that of the first case. 

In Fig. 1151 we show the lowest-lying 16 eigenvalues on 10 gauge configurations for each of the above two cases. Every 
line in these plots connects all the 16 eigenvalues in one configuration. It is evident that on the finer lattice, the 
4-fold degeneracy appears, while it is much less clear for the coarser lattice. Although taste-symmetry restoration 
appears on fine lattices in our work, from Fig. 1151 it is indicated that such restoration does not show up on the 
coarse lattices in our simulations. Such effects are expected, since unimproved staggered fermions are implemented. 
This necessitates good control of the continuum extrapolation, which is addressed in detail in Sec. IV Dl 
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0.7810709(57) 


6 


11.15 


0.8160114(69) 


6 


13.85 
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20 
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20 
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20 
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20 
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20 
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20 
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20 
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20 
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20 


50.00 


0.9597892(11) 



TABLE III: The expectation values of the plaquette in a significant fraction of our simulations. 



Appendix B: Values of plaquette and TPL coupling constant raw data 



In this appendix, we present details for the plaquette values and the TPL scheme renormalised coupling constants 
obtained at various lattice volumes, L/a, and bare couplings, /3. We also give the numbers of HMC trajectories for 
the computation of the TPL coupling. 
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TABLE IV: Raw data for the renormalised coupling in the TPL scheme. 
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